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We provide a detailed derivation of the low-energy model for Zn-diluted La2Cu04 in the limit 
of low doping together with a study of the ground-state properties of that model. We consider 
Zn-doped La2Cu04 within a framework of the three-band Hubbard model, which closely describes 
high-T c cuprates on the energy scale of the most relevant atomic orbitals. To obtain the low-energy 
effective model, we first determine hybridized electronic states of CuCU and Zn04 plaquets within 
the Cu02 planes. Qualitatively, we find that the hybridization of zinc and oxygen orbitals can 
result in an impurity state with the energy e, which is lower than the effective Hubbard gap U. In 
the limit of the effective hopping integral t <e, U, the low-energy, spin-only Hamiltonian includes 
terms of the order t 2 /U and t 4 /e 3 . That is, besides the usual nearest-neighbor superexchange J- 
terms of order t 2 /£/, the low-energy model contains impurity-mediated, further-neighbor frustrating 
interactions among the Cu spins surrounding Zn-sites in an otherwise unfrustrated antiferromagnetic 
background. These terms, denoted as J% n and J z ; n , are of order t 4 /e 3 and can be substantial when 
e ~ U/2, the latter value corresponding to the realistic CuCb parameters. In order to verify this 
spin-only model, we subsequently apply the T-matrix approach to study the effect of impurities on 
the antiferromagnetic order parameter. Previous theoretical studies, which include only the dilution 
effect of impurities, show a large discrepancy with experimental data in the doping dependence of 
the staggered magnetization at low doping. We demonstrate that this discrepancy is eliminated by 
including impurity-induced frustrations into the effective spin model with realistic CuCb parameters. 
Recent experimental study shows a significantly stronger suppression of spin stiffness in the case of 
Zn-doped La2Cu04 compared to the Mg-doped case and thus gives a strong support to our theory. 
We argue that the proposed impurity-induced frustrations should be important in other strongly 
correlated oxides and charge-transfer insulators. 

PACS numbers: 75.10.Jm, 75.40.Gb, 75.30.Ds, 75.50.Ee 



I. INTRODUCTION 

Impurity effects in spin systems have attracted consid- 
erable attention since the 1960s as a part of the broader 
interest in defects in solids^E' and also in the context of 
percolative phenomenaP^to which diluted magnets pro- 
vide excellent experimental realizations. More recently, 
the discovery of high-T c superconductivity in insulating 
antiferromagnets doped by mobile charge carriers has 
lead to an extensive research effort in various aspects 
of the problem, including role of d isorde r in quantum 
critical states,^ quantum percolation* 10 1 11 1 and others. It 
has also been realized that impurities may serve as valu- 
able local probes into the properties of strongly correlated 
systems in general, revealing impo rtant aspects of their 
electronic degrees of freedom 

Out of the parent compounds of cuprate supercon- 
ductors, it is La2Cu04 that has been studied most 
comprehensivelyP^HS] In its pristine form, it is an excel- 
lent realization of the two-dimensional, spin-^, nearest- 
neighbor square-lattice Heisenberg antiferromagnet 22 " 24 * 
formed by Cu 2+ ions surrounded by oxygens. The dilu- 
tion is achieved by chemically substituting isovalent Zn 2+ 
or Mg 2+ spinless ions for Cu 2+ , see Fig.yja). Thus, it is 
only natural to expect that the proper low-energy model 
of La2Cui_ :c (Zn,Mg) :c 04 must be the nearest-neighbor 
site-diluted Heisenberg mo del see Fig.JlJb). In order 
to elucidate the properties of La2Cu04 diluted by spin- 
less impurities, and of the associated site-diluted spin 




FIG. 1. (Color online) (a) A schematic view of Zn-doped 
Cu02 plane, (b) Site-diluted Heisenberg model representa- 
tion of it. The arrows are Cu spins and the lines denote su- 
perexchange interactions, (c) The same, with extra frustrat- 
ing interactions between the next- (Jzn) an d the next-next- 
(Jzn) nearest neighbor Cu sites, surrounding Zn impurity. 



model, extensive experimental studies have been per- 
formed using nuclear magnetic resonance (NMR) [nu- 
clear quadrupole resonance (NQR)], muon spin relax- 
ation (/iSR), elastic and inelastic neutron scattering, and 
magnetometry.^^U An equally comprehensive theoret- 
ical effort included the spin-wave T-matrix, Quantum 
Monte Carlo, and numerical real-space 1/S calculations 
of a number of quantities that allowed for extensive cross- 
examinat ions pSHSHl 

One quantity in particular, the average magnetic mo- 
ment per Cu, M(x), has been investigated in detail. The 
doping dependence of M(x) is the purely quantum effect 
related to the impurity-induced suppression of the or- 



2 



der parameter because the normalization to the number 
of magnetic ions naturally separates the classical effect 
of dilution from it. While the overall results for M vs 
x show a reasonable agreement, a substantial discrep- 
ancy between the experiment and the theory, both ana- 
lytical and QMC, has been observed. The experimental 
results indicate a substantially stronger — a factor of ap- 
proximately 2 — suppression of the order parameter M 
per impurity due to disorder-induced quantum fluctua- 
tions, and the experimental data for M{x) are also always 
below the theoretical curves pSEU 

In our recent work, 39 a resolution to this problem has 
been suggested: impurities should not be considered as 
electronically inert vacancies that simply eliminate inter- 
actions among surrounding Cu spins. In addition to the 
dilution, the hybridized electronic states of the impurity 
and of the nearest oxygens can provide extra degrees of 
freedom that generate longer-range frustrating interac- 
tions, schematically shown in Fig. [TJc). Such impurity- 
induced frustrating interactions can be expected to sig- 
nificantly enhance local quantum fluctuations. In Ref.l39l 
we have outlined our approach to the problem and pre- 
sented our results for M(x) together with the comple- 
mentary QMC results. The latter have unequivocally 
supported the same conclusion: the dilution-frustration 
model exhibits stronger suppression of the order and, for a 
choice of parameters appropriate for Zn-doped La2Cu04, 
bridges the gap between experiments and previous the- 
oretical calculations based on the dilution- only model. 
Since then our theory has received further experimen- 
tal confirmation from the /iSR studies of spin stiffness 
in Zn and Mg doped I^CuCUp^ which demonstrated a 
stronger suppression of spin stiffness in Zn-doped case, 
in agreement with the expectations from the theory. 

In this work, we expose the details of the derivations 
of the dilution-frustration model and of the subsequent 
analytical calculations of the doping-dependent magne- 
tization. The key element of our theory is that it starts 
from a realistic model, the site-diluted three-band Hub- 
bard model, instead of already "effective" models of the 
Cu02 plane, such as the diluted Heisenberg or the diluted 
one-band Hubbard model. 

We begin with Zn-doped I^CuC^ within the frame- 
work of the impurity-doped three-band Hubbard model, 
which closely describes the diluted high-T c cuprates and 
other transition- metal oxides on th e energy scale of the 
most relevant atomic orbitals l- 41 * 42 * Similar to the deriva- 
tion of the one-b and H ubbard model from the three-band 
Hubbard model ] 43 * 44 ! the cell-perturbation approach is 
used to describe hybridization of the energy levels of Cu 
and Zn with oxygen orbitals. The approach does not 
require Cu-0 and Zn-0 hopping integrals to be smaller 
than the charge-transfer gap, therefore the locally hy- 
bridized states on Cu or Zn and surrounding O's are di- 
agonalized without approximation!^" 4 ^ Then the three- 
band model can be rewritten as a "multi-orbital" Hub- 
bard model with the effective "Cu" and "Zn" states con- 
nected by effective hoppings. Since the structure of the 



lowest states in the single-particle and two-particle sec- 
tors of the multi-orbital model is the same as in the one- 
band Hubbard model {i.e. the lowest two- hole state is 
the Zhang- Rice-like singlet), the equivalence of the two 
models can be justified.^ED i n this approach, for the 
dilution-only picture to be valid the effective "Zn" -states 
must not occur below the effective Hubbard J7, the situa- 
tion more likely to be valid for the case of Mg-doping due 
to the lack of available electronic states on Mg ion. For 
the case of Zn-doping, its electronic states^BEU hybridize 
with the states on the oxygen orbitals and can result in 
the states below the Hubbard energy gap. This provides 
surrounding Cu-spins with additional virtual states to 
execute their superexchange processes through, thus fa- 
cilitating extra couplings that connect spins in the imme- 
diate vicinity of impurity. Therefore, the spinless impu- 
rity, in effect, can lead to a cage of frustrating interactions 
around itself, shown in Fig. [TJc). 

For the sake of a qualitative picture, and in the spirit 
of mapping of the multi-band Hubbard model onto the 
single-band one, the result of this consideration is that 
the impurity-doped system is not equivalent to the site- 
diluted Hubbard model with electronically inert impurity 
sites, but rather to the t-e-U model, where in addition to 
the usual hopping t and the two-particle energy gap U 
there is the lowest energy state of the effective impurity 
sector, denoted as e. At half- filling, the t-U part reduces 
to the Heisenberg model with J~t 2 /U, as usual, with the 
higher-order terms ~t A /U 3 negligible if t <C Uf^ieM Vir- 
tual transitions through the effective impurity level e gen- 
erate superexchange interactions of the order of t 4 /e 3 be- 
tween the next- and the next-next-nearest-neighbor Cu- 
spins ( Jz n and Jz n ) that are also nearest-neighbors of the 
impurity site. When the energy at the impurity site e is 
less than the Hubbard gap, such terms are not negligible 
and may be comparable to J. For an estimate, taking 
e = U/2 and U/t = 10 gives J' Z J J ~ (t/U) 2 (U/e) 3 -0.1, 
and, given the geometry of the square lattice, the com- 
bined effect per impurity is J|°* = 4J^ n +2 J£ n -0.6 J. The 
origin of J% n and J% n is clearly distinct from the gener- 
ally considered next- and next-next-nearest-neighbor su- 
perexchange interaction (J2 and J3), which are of the 
order of t A /U 3 and do not affect the order parameter 
specifically due to dilution! 5 - 5 ^ 

While the technical details and the effort of the present 
work are more involved, this qualitative t-e-U model 
properly reflects the key idea of our approach. In 
practice, we perform a similar type of the 4 t/l -order 
expansion 54 of the multi-orbital Hubbard model, keep- 
ing track of all the relevant one- and two-hole states 
of the model and their dependencies on the original 
three-band model parameters. This allows us to per- 
form a detailed microscopic calculations of J% n and J![ n 
and estimate their values. Since the experimental value 
of the nearest-neighbor superexchange for La2Cu04 is 
known (J~0.13eV), it can be used to narrow down the 
range of parameters of the three-band model as was done 
previously! 4 ^ Although the electronic parameters of Zn 
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states, such as the energy of the bare Zn- level and the 
Zn-0 hybridization, are not known precisely ] 49 * 5 1 * we vary 
them to verify that the energy of the lowest impurity level 
e indeed falls comfortably below the effective Hubbard 
U for a wide and reasonable range of both parameters. 
By projecting the multi-orbital Hubbard model onto the 
low-energy spin-only model, we analyze possible value of 
the total frustrating effect per impurity and find it to 
be between J|°* ~ (0.2 — 1.0) J for the same range of 
parameters. 

In addition, we also provide a detailed analysis of 
the individual processes that contribute to J' Zn and J% n . 
Specifically, we have found that, counterintuitively, the 
longer-ranged J% n is greater than J% n . This is due to 
a subtle cancellation between the "regular" 4 t/l -order 
superexchange processes and the analogs of the ring- 
exchange-type processes involving Zn and three Cu sites 
on a nearest-neighbor plaquette, see Fig. [TJa), with the 
latter contributing to J^, but not to J% n . That results 
in a stronger frustrating coupling between copper spins 
across the Zn-site, with the ratio J% n / 'J% n c± 2 ~ 4 in a 
wide range of the three-band model parameters. Alto- 
gether, the calculation based on the three-band model, 
albeit more involved, provides a strong support to our 
central idea and gives an order-of-magnitude estimate of 
the frustrating terms in the dilution-frustration model. 

Upon establishing the structure of the low-energy spin- 
only model for the diluted system, to which we refer to 
as to the dilution-frustration model, we investigate the 
impact of impurities in this model on the order param- 
eter M . This is achieved by means of the analytical 
T-matrix approach within the spin-wave approximation 
based on the exact diagrammatic treatment of the impu- 
rity scattering amplitudes and subsequent disorder av- 
eraging. Technically, we closely follow the approach of 
Refs.|33|and|3l 

First, we apply the spin- wave approximation to rewrite 
the dilution-frustration model on the 2D square lattice. 
After the Fourier and Bogolyubov transformation we de- 
compose the impurity scattering matrices into the s-, 
p-, d-wave orthogonal components with respect to the 
scattering site. Then, the T-matrix approach is used to 
solve exactly the problem of scattering off one impurity in 
each of the scattering channels. The subsequent disorder- 
averaging approximates impurities as independent ran- 
dom scatterers and effectively restores translational in- 
variance for spin excitations propagating in an effec- 
tive medium. Such an approximation neglects impurity- 
impurity interaction effects, which are expected to be 
small at small doping. Disorder averaging extends the 
T-matrix approach to the finite impurity concentration 
and yields the spin-wave self-energies as simply related 
to the forward scattering components of the T-matrix. 
Next, for the given impurity concentration and values 
of frustrating parameters, the on-site ordered magnetic 
moment M is calculated from the renormalized magnon 
Green's functions. 

Compared to the dilution-only model, the modification 



of this method for the dilution-frustration model con- 
cerns changes in the p- and <i-wave scattering channels, 
while the s-wave channel can be shown to be unaffected 
by the frustrating terms. Generally, the advantage of 
the T-matrix method is that it offers a systematic way 
of studying the order parameter M as a function of the 
concentration x and of the parameters J% n and J% n . 

One of the unusual findings in this study is that the 
next-next-nearest neighbor J£ n frustrating bond sup- 
presses the order as effectively as two next-nearest J' Zn 
bonds of the same strength. This result has also been 
supported by the QMC calculations, as was discussed 
previously.^ Given that according to the three-band 
model calculations the J^-term is larger than the J% n - 
term, this finding underscores its importance for the 
mechanism of impurity-enhanced suppression of the or- 
der parameter. 

We find that the experimental rate of suppression of 
the order in Zn-doped La2Cu04 is met by the T-matrix 
results of the dilution- frustration model at J|°* = 0.28 J 
if we fix the ratio of the two frustrating terms to 2, 
J'L = 2J L = 0-07J, or at J|° n * = 0.24J for J£ n = 
4Jz n = 0.08J, according to the three-band model results. 
As was discussed in our previous work, 39 QMC results 
seem to suggest a somewhat higher value ^ 0.4 J 
(Jz n = 27^ > 0.1 J), although they are obtained from 
the finite-size extrapolations that may overestimate J%£. 
Both the T-matrix and the QMC results correspond to 
a modest amount of frustration, well within the window 
suggested by the three-band model calculations. With 
our analytical and numerical results agreeing quantita- 
tively with each other, we have suggested further high- 
precision experiments at low doping. One of such exper- 
imental confirmations has come recently from the /iSR 
studies of Zn- and Mg-doped La2Cu04,^ which has 
shown a substantially stronger suppression of the spin 
stiffness in the case of Zn-doping, in agreement with the 
expectations from our theory. 

The paper is organized as follows. In Sec. [H] we discuss 
the diluted three-band Hubbard model within the cell- 
perturbation approach and derive the multi-orbital Hub- 
bard model from it. We study the hybridized impurity 
level and analyze its energy with respect to the effective 
Hubbard gap for a range of three-band model parame- 
ters. We derive an effective low-energy spin-only model 
by applying canonical transformation to the multi-orbital 
Hubbard model at half-filling. The virtual steps leading 
to the frustrating superexchange interactions, J% n and 
Jz n , are analyzed and their values are calculated. In 
Sec. Ill, we study the effective dilution-frustration model 
using the T-matrix approach and disorder averaging. We 
derive the staggered magnetization as a function of im- 
purity concentration and frustrating couplings. The re- 
sults are compared with experimental data. Section IV 
contains our conclusions. Appendices A and B contain 
details of the three-band model and the T-matrix calcu- 
lations, respectively. 
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II. EFFECTIVE MODEL OF COPPER-OXIDE 
PLANE WITH ZINC IMPURITIES 

In this section, we derive an effective low-energy model 
for the Zn-doped La2Cu04 at half-filling. We be- 
gin with the consideration of the realistic, Zn-diluted 
three-band Hubbard model, which contains additional 
impurity states associated with Zn. Using Wannier- 
orthogonalization for O-orbitals and a more natural lan- 
guage of the locally hybridized CUO4 and Zn04 states, 
we first rewrite the three-band Hamiltonian as the multi- 
orbital Hubbard model. We discuss the new feature of 
the model — the hybridized Zn and O orbitals can result 
in an impurity state with the energy that is lower than 
the effective Hubbard gap. The subsequent transforma- 
tion to the low-energy, spin-only model is known as the 
cell-perturbation method. The local basis of CUO4 and 
Zn04 states of the multi-orbital Hamiltonian provides a 
natural small parameter for such a projection, the effec- 
tive hopping between the states of the nearest-neighbor 
clusters. We proceed with this transformation and give 
a quantitative analysis of the individual processes that 
contribute to the superexchange terms of the effective 
model. Next, we analyze possible range of the frustrat- 
ing interactions in the low-energy model of the Zn-doped 
Cu02 plane. The discussed approach should be valid as 
long as the system is in the Mott insulating state. 



A. Three-band to multi-orbital Hubbard model 

The three-band Hubbard model, which was proposed 
to describe relevant electronic degrees of freedom of the 
Cu02 planes of the high-T c cupratesj^ is formulated in 
terms of the hole creation and annihilation operators act- 
ing on the vacuum state of the completely filled 3d 10 
shells of Cu and 2p 6 shells of O 



U = € d U la + £ P Yl U 
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where e& and e p are the energies of the hole in the copper 
d x 2_ y 2 and in the oxygen p x or p y orbitals, respectively. 
The copper (oxygen) sites are labeled with the index Z(m), 
the number operators are nf a = d] a d la (n p ma = pLjw)> 
and a =t, 4- is the spin. The Hubbard repulsion in the 
copper d x 2_ y 2 orbitals is Ud and the hopping between 
copper and oxygen orbitals is t p d- For the Cu-0 hoppings 
we use the convention 42 in which the relative signs of or- 
bitals are absorbed into the definition of p ma (p ma ) oper- 
ators. The summation (Im) is over the nearest-neighbor 
Cu-0 bonds. 

The relevant transitions within the three-band model 
involve d 10 -d 8 states on copper and p states on oxygen 
sites. At half-filling, there is one hole per Cu02 unit cell 
and the ground state is the antiferromagnetic insulator.^ 
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FIG. 2. (Color online) Zn-doped CuCb plane with natural 
partitioning in Cu04 and Zn04 clusters. 



The localized 5 = 1/2 spins, forming the Neel-ordered 
state on the square lattice, are provided by the holes 
that are predominantly in the d 9 states on Cu and are 
hybridized with the p 5 states on O, see Fig. [2] 

The minimal form of the three-band model in ([!]) is 
often supplemented with additional terms, such as the 
direct oxygen-oxygen hopping, on-site Coulomb interac- 
tions on O-sites, and the nearest-neighbor repulsion be- 
tween O- and Cu-holes. The main effect of such terms is 
in a quantitative renormalization of the energy levels and 
hopping integrals) 44 1 46 1 leaving the results obtained with 
([I]) qualitatively the same. 

The substitution of the isovalent Zn 2+ ion with the 
nominally completely filled 3<i-shell for Cu 2+ leaves the 
oxygen lattice translationally invariant, see Fig. [2] While 
the electronic levels of Zn doped into a Cu02 plane have 
not been determined preciselyJ^HSU it is generally agreed 
that the relevant states may occur in a reasonable vicinity 
of the oxygen level, Sz n — %> — 2 — 5eV, and we treat this 
energy as an adjustable parameter in the following. Thus, 
the impurity Hamiltonian can be written as 



5H Zn = e Zn < + U Za £ nffnff 

ta I 

UKX 

+ H.c), (2) 

(£m)a 

where Zn sites are labeled with £, szn is the energy of 
the relevant orbital on Zn, Uz n is an effective Hubbard 
repulsion in that orbital, and tznO is the hopping integral 
between nearest-neighbor Zn and O sites, yet another ad- 
justable parameter. The hole creation/annihilation op- 
erators on Zn are a\ a and a ia and nf£ — a \ a a £a' 

The checkerboard structural motif of the Cu02 plane 
in Fig. [2] suggests a natural partitioning of the localized 
states into symmetric CUO4 and ZnC^ clusters, in which 
each Cu or Zn is surrounded by four oxygens The ad- 
vantage of such basic units is in constructing the local 
states that allow to take into account local hybridiza- 
tion of Cu or Zn and surrounding O states without ap- 
proximation, while the remaining hybridization between 
the clusters can be treated perturbatively. In particular, 
such a cell-perturbation approach does not require Cu-0 
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and Zn-0 hopping integrals to be much smaller than the 
charge-transfer gap A = e p — s^P™^ 

On the other hand, the linear combinations of the oxy- 
gen states in CUO4 (Zn04) cluster are not orthogonal 
to the ones in the nearest-neighbor clusters. The ele- 
gant Wannier-orthogonalization procedure, suggested in 
Ref. [43j treats the O-lattice separate from the Cu and 
Zn, and, via orthogonalizing p-operators in the k-space, 
leads to the basis of the symmetric-oxygen orbitals 



{P x m 



(3) 



that are now associated with the same site index as the 
Cu (Zn) of the CUO4 (ZnC^) cluster. The details of this 
procedure are given in Appendix [Aj 

With that, it is natural to divide the Cu-0 (Zn-O) 
hopping terms in and Q into the local part, which 
corresponds to hybridization within one cluster, and into 
the hopping part, corresponding to the coupling between 
the states in different clusters. We note, that due to 
the non-local nature of the Wannier-orthogonalization, 
the hopping parts now contain terms that go beyond the 
nearest-neighbor clusters. Thus, the Hamiltonian of 0, 
written using the orthogonalized basis of O-states (|3| is 



H loc = £ L d nf a + e p nl + %f a nf a 

la 1 1 

-2X t pd (4 a q la + K.c.) J 



(4) 
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H.c. 



wher e A/// are the Wannier amplitudes given in 
(A7). First, the amplitudes \n> decrease rapidly with 
distance,^ \n> ~ l/|r/ — i> | 3 , so the terms beyond the 
nearest neighbor in the hopping part of Q can be ne- 
glected. Second, the intra-cluster amplitude Ao = 0.9581 
is close to 1, and thus is taking most of the hybridization 
into account. The nearest-neighbor amplitude \(w) = 
0.1401 leads to the effective hopping parameter between 
different states that is significantly reduced compared to 
the bare t p d hopping. This provides a strong justification 
for the subsequent perturbative expansion in such an ef- 
fective hopping relative to the effective Hubbard gap, the 
latter largely defined by the local hybridization. 

Applying the same transformation of the O-states ([3| 
to the Hamiltonian of the Zn-cluster ([2| we obtain 

lot ^ 

-2A t Z no(aL^ + H - c -)} (5) 
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Note that even if Zn-states are completely neglected 
in ffl^n? Zn04 cluster still has an unoccupied oxygen 
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FIG. 3. (Color online) (a) A schematic view of the energy 
levels Ei and Ei of one- and two-particle sectors in Cu and 
Zn clusters in Eqs. (|6| and |7]). (b) The levels of the effective 
t — e — U Hamiltonian that keeps the lowest energy states in 
each sector. 



state^ with the energy e p , which is higher than the hy- 
bridized CUO4 states of the surrounding clusters, and it 
permits virtual superexchange processes through it. 

The next step in the cell-perturbation approach is the 
diagonalization of H loc and 57-L loc . Such a diagonaliza- 
tion is performed separately for the states with different 
number of holes in a cluster. Because the system is close 
to the half-filling, the most relevant states are the one- 
and two-hole states, and the states with more holes are 
higher in energy, see Fig. [3ja). The full set of di, an, and 
qi one- and two-hole states of Cu- and Zn-clusters are 
listed in Appendix [A| (|A8l)-( |A13[ ), see also Ref. HU 

After the diagonalization of the local parts of the 
Hamiltonian, the hopping terms are also rewritten in the 
basis of the new states. Altogether, the three-band model 
H of @ can be rewritten as a "multi-orbital" Hubbard 
model with the effective "Cu" eigenstates with local en- 
ergies Ei connected by effective hoppings F-3 : 



+ EE^' (m\4)(4\(t 

(W) ii'jj' 

Similarly, for 5H of ^ 

5u = Y,m\)m 



(6) 



H.c.) . 



(7) 



EE^' (i^)i^')w"k^i+h.c.) . 
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The hybridized, orthogonal sets of local states are |^) 
for CUO4 and for Zn04 clusters, where l(t) is the 
site-index and i is labeling the states. For example, the 
three one-hole states in Cu cluster in Fig. [3ja) are the 
"bonding" and "antibondi ng" mix o f Cu and symmetric 
O states d\ and qi : Eqs. ( |A4[) and (A8), and the anti- 
symmetric O-state (qi, Eq. ( |A4[ )). The hopping integrals 

between adjacent clusters, and F?P, , connect initial 
and final states z, i r and j, j f '. 



The Hamiltonian in |6| and ^ is "multi-orbital" be- 
cause each n-particle sector contains more than one state. 
However, since such states are locally orthogonal and we 
are interested in the low-energy effective model of (|6| and 
([7]), the most important transitions involve the lowest 
states from each sector, see Fig. |3ja). Importantly, the 
structure of the lowest states in the single-particle and 
two-particle sectors of the multi-orbital model is the same 
as in the one-band Hubbard model, justifying the close 
correspondence of the three-band and one-band Hubbard 
modelsP3E3 The lowest one-hole states in Cu- and Zn- 
cluster are S = 1/2 doublets, the linear combinations 
of the oxygen and the Cu(Zn)-orbitals, while the low- 
est two- hole state is the Zhang- Rice-like singlet, a mix 
of th r ee di fferent singlets [Cu-Cu, 0-0 and Cu-O], see 
|A8|)-(|AT2|). 



The new physics is brought in by the possibility of 
an unoccupied impurity level, the lowest state from the 
"Zn" single-particle sector, to be lower than the effective 
Hubbard U on the "Cu" sites. This leads to an effective 
t — e — U model, see Fig. [3^b), which provides surround- 
ing Cu-spins with additional virtual channel for the su- 
perexchange processes, connecting spins in the vicinity 
of impurity. Conversely, if the effective "Zn" -states do 
not occur below the Hubbard gap (the case of Mg 2+ due 
to the lack of available states on it), the impurity can be 
treated as electronically inert, leading only to dilution, 
but not to extra interactions. Parameters that control 
the impurity level are discussed next. 




FIG. 4. (Color online) (a) Energy of the lowest state from the 
Zn-cluster single-particle sector, vs the energy difference 
between the bare one-hole levels on Zn and O, Az n = £zn-^p, 
for two values of tznO- (b) £gff vs tzno/t p d for two represen- 
tative values of Az n . Parameters of the three-band model are 
fixed as described in text, t p d = 1.5eV, A = 3eV, £/d^>A. The 
resultant Hubbard gap is U e s ~ 3.6eV. 



values of tznO? and in Fig. Qb) on the the hybridization 
tznO f° r two representative values of Az n . We also show 
the effective Hubbard energy /7 eff to demonstrate the va- 
lidity of the qualitative level structure in Fig. [3] for a wide 
range of parameters. Altogether, the electronic levels of 
Zn and their hybridization with O-states are important 
in lowering e^§, and is below the Hubbard gap for 
the realistic parameters of the model. 



B. Parameters and effective impurity level 

The realistic values of the three-band model parame- 
ters for the Cu02 plane have a certain degree of vari- 
ation because none of them is measurable directly and 
they come as a result of a parametrization of the first- 
principles calculations. Thus, the Hubbard repulsion on 
Cu is Ud = S— 12eV, the charge-transfer gap A = S p— Sd = 
2-4eV, and the Cu-0 hopping is t pd = 1 - 1.5eV P ' 4 ^ ' 47 ' 
A physical approach to fix them to a particular set of 
values was suggested in Ref. [471 One can require that 
the observables, such as the Cu-Cu superexchange or the 
optical gap, calculated from the model (|6| match their 
observed values. In our case, we fix them to Ud(Uz n ) ^ A, 
A = 3eV, and t p d = 1.5eV to yield the experimental 
value of the superexchange J ~ 0.13eV in the effective 
low-energy model discussued in Sec. [H|C. Among other 
things, the effective Hubbard gap for the Zhang-Rice- 
like singlet and the ground state half-filled system is 
Z7 eff ^3.6eV. 

The parameters of Zn states, such as the bare one- hole 
level on Zn, Az n = £z n -e j9 , and the Zn-0 hopping tznO? 
are not precisely knownPSHHI As is shown in Fig. [ZJ we 
vary them substantially to determine the resulting energy 
of the lowest level in the single-particle sector of the Zn04 
cluster, e^Q. In Fig. [4^a) we show its dependence on the 
energy of the bare Zn-level, Az n , for two representative 



C. Projecting to the spin-only model 

After establishing the validity of the t-s-U -like level 
structure of Fig. [3] for the realistic parameters of the 
three-band Hubbard model with Zn impurity, we now 
turn to the low-energy properties of the model ([6| and 
([7]). To derive the low-energy model one needs to con- 
sider virtual transitions between different states. As was 
argued above, since the system is at half- filling, the rele- 
vant transitions are between the lowest states from one- 




Cu Cu Zn Cu 



FIG. 5. (Color online) Schematic view of various virtual hop- 
ping processes between relevant energy levels of the model Q 
and (m), see text for notations. 
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FIG. 6. (Color online) Superexchange processes for (a) 5Ji, 
(b) 5J2. See text and Fig. [5]for notations. 
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FIG. 7. (Color online) Ring- exchanges for (a) 6J3, (b) SJ4. 



and two-hole sector of both Cu- and Zn-cluster states. 
Altogether, there are five different hopping integrals, as 
is shown in Fig. [5] For brevity, we switch to £'s for the 
hopping integrals from and in and |?|. For 
instance, t\ = F^q 1 is the hopping of the spin to the neigh- 
boring empty site and £2 = F?i is creating the Zhang- 
Rice singlet from the two one-hole sites leaving the other 
site empty. For the transitions involving Zn- states, tn 
is between one-hole states on Cu and Zn, while ti2 and 
£21 are between the one- hole and two- hole states on Cu 
and Zn. The higher-energy transitions will be neglected 
as their contributions are small. The explicit expressions 
for hopping integrals can be obtained by evaluating ma- 
trix elements of the Hamiltonian in terms of the original 
Cu, Zn, and O operators, Q and ([5]), between the initial 

and final states in the basis of the local eigenstates 
and see Appendix |a| for details. 

Next, we apply a canonical transformation to this ex- 
tended t-e-U-like model assuming that the hopping in- 
tegrals are smaller than the effective and s^. At 
half- filling and in the second order, the t-U part yields 
the Heisenberg model with J = 4|t 2 | 2 /^ff U and with the 
negligible higher-order terms. In the presence of impu- 
rity, the same transformation leads to the dilution-only 
model, in which the four adjacent spin-spin links are cut 
by the spinless Zn-site. For the sites in the vicinity of 
impurities, we extend the transformation to the fourth 
ordei^M^in t/e(U) to include the -e- part of the model 
for the purpose of taking into account the effects of the 
in-gap impurity state. 

In addition to deriving the low-energy, spin-only model 
with the parameters that can be traced to the three-band 
ones, we are also able to analyze individual superex- 
change processes that contribute to the terms of that 
model. We find that in the 4 th order there are two types 
of virtual transitions through the impurity site. First 
involves three sites, two Cu and one Zn, with the cop- 



pers either across the Zn-site or at the right angle, as in 
Fig. [6j Second type needs three Cu in addition to Zn 
impurity, arranged in a 4-site plaquette, see Fig. [7] Each 
of the transition types yields two superexchange chan- 
nels shown in Figs. [6] and [7| First two, in Fig. |6ja) and 
(b), are the standard superexchanges, taking four vir- 
tual steps instead of the usual two. We denote the cor- 
responding couplings generated by such processes as 8J\ 
(doubly-occupied state on Cu) and 6 J2 (doubly-occupied 
state on Zn), respectively. The other two processes, in 
Fig.^a) and (b), are the ring exchanges, which we de- 
note as 5 Js (no doubly-occupied state involved) and SJ4 
(doubly-occupied state on Cu), respectively. These cou- 
plings are given by 

SJl = 4| ^ 2l|2 2 > *J* = 8Ml2 ' 2 2 , (8) 

U^ielif ([/^ + 2 £ zn)( £ zn) 2 

XJ 2|t 11 t 1 | 2 4|tii*i* 2 *2l| 

in terms of the notations of Fig. [5] Assuming that e 7 ^ is 
smaller than Z7 e ff according to the discussion of Sec.[ll|B, 
the magnitudes of SJ^s should substantially exceed the 
conventional fourth-order terms of order ^ff/^eff' which 
we neglect. Another important distinction of the cou- 
plings in (J8|, is that they occur due to impurities and thus 
must have a direct relation to the doping-dependent ef- 
fects. Altogether, the considered processes combine into 
the next-nearest- and next-next-nearest-neighbor frus- 
trating interactions around Zn-impurity, see Fig. [TJc), 

JL = SJ 1 +SJ 2 + SJ S + SJ 4 , (9) 

It is clear that only superexchange-type processes of 
Fig. [6] can contribute to the next-next-nearest-neighbor 
couplings across the impurity J^, while the ring- 
exchanges of Fig. [7] are also at play for J' Zn . Importantly, 
the latter terms have a ferromagnetic sign, reducing the 
value of J% n compared to J% n . Thus, counterintuitively, 
the longer-ranged J% n should be greater than J' Zn . This 
feature is made explicit in our Fig. |8j which shows that 
the ratio of J% n / / J^ n varies from about 4 to 2 for a rep- 
resentative range of the electronic parameters of Zn with 
the other three-band model parameters as in Fig. [4] 
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FIG. 8. (Color online) The ratio of Jzn/^Zn vs tz n o/t p d for 
two representative Az n - Other parameters are as in Fig. [4] 
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FIG. 9. (Color online) The total frustration effect per Zn 
impurity = 4J Zn + 2 Jz n (a) vs Az n for several represen- 
tative values of tznO and (b) vs tzno/tpd for several values 
of Azn, respectively. Three-band parameters are £ p d = 1.5eV, 
A = 3eV, Ud^> A as before. The shaded area is the range 
of J Z n* needed to explain experimental re duct ion of the stag- 
gered magnetization M(x)/M(0), see Sec. Ill 



As is shown in Fig. [TJc), each impurity generates four 
Thus, the total frustrating effect per 
2Jz n . We show the depen- 



and two J Zr] 



Zn-impurity is J^f 
dence of J*g 



9lwith the same set 
TUB and in Fig. HI 



on A Zn and t Zn o in Fig. 
of three-band parameters as in Sec. 
The shaded area in the graphs shows the range of 
needed to explain experimentally observed reduction of 
the staggered magnetization M(x)/M(0), according to 
the T-matrix discussion of Sec. |III| Clearly, the uncer- 
tainty in the electronic parameters of Zn-orbitals does not 
allow us to provide significant restrictions on the value 
of J|n £ - However, the values needed for M(x)/M(0) dop- 
ing dependence outline the minimal requirements on such 
three-band model parameters. In particular, the cou- 
pling of Zn and O orbitals should not be much weaker 
than that of Cu and O, while the position of the level 
on Zn is much less restricted. Overall, the realistic re- 
quirements on Zn electronic degrees of freedom leave a 
wide range of J|°* and they support the validity of our 
proposed impurity-induced frustration mechanism. 

Finally, the effective low-energy, spin-only model for 
I^Cui-^Zn^C^ can be divided into two parts 

T~L = J ^2 ^ l ' ^ ^zn Si ■ Sit + Jz n ^2 Si • Sj/, 

(U') (n% (ii'Yl 

(10) 

the dilution-only model (first term), and the impurity- 
induced frustrating terms, with and (W)" denoting 
next- and next-next-nearest-neighbor bonds that are also 
nearest-neighbors of the impurity site i, see Fig. []Jc). 
The summations are carried over the Cu sites only. 

We propose that the dilution-frustration model ( |T()| ) 
provides a proper description of Zn-doped Cu02 planes 
and discuss its properties next. 



III. ON-SITE MAGNETIZATION IN THE 
LOW-ENERGY MODEL 

With the structure of the low-energy spin-only model 
for Cu02 plane diluted by Zn impurities given in ( 10), we 
investigate in detail the impact of impurities on the or- 
der parameter, the average on-site magnetization M. We 
use the linear spin- wave approximation to the problem of 
single impurity in the model (10) and apply an analytical 



T-matrix approach to it, which is based on summing up 
exactly infinite diagrammatic series for the impurity scat- 
tering amplitudes. Then we approximate the problem of 
finite concentration of impurities as a linear superposi- 
tion of scattering effects of individual random impuri- 
ties, which amounts to an "effective medium" approach 
via disorder averaging. This step restores translational 
invariance for magnons but modifies their dispersion and 
provides them with damping through the impurity scat- 
tering self-energies. The averaging also takes into ac- 
count the increase of the population of magnons around 
impurities, i.e., for the extra fluctuations of spins that 
reduce the value of the ordered moment. The latter ef- 
fect is calculated from the renormalized magnon G reen' s 
function as a function of impurity concentration x! 33 * 34 * 

We note that the order parameter is normalized by the 
number of magnetic ions: M(x) = J2i l^l/^m, where 
N m = N — Ni mp is the number of magnetic sites (Cu 2+ ) 
and N is the total number of sites. The reason for that 
is twofold. First, the experimental data on M(x) from 
such techniques as NMR[NQR] and neutron scattering 
are naturally normalized this way. Second, such a nor- 
malization also separates the classical effect of dilution 
from the purely quantum-mechanical suppression of the 
order parameter. For instance, M(x) defined this way 
is close to a constant for the classical (Ising) antiferro- 
magnet on a square lattice, as is shown in Fig. |l0|(a), 
because this quantity is equivalent to a probability for a 
spin to belong to an infinite cluster, which is very close 
to 1 below the percolation threshold, x p w 41%. 

In this Section, we first demonstrate the discrepancy of 
the previous theoretical works based on the dilution-only 
model with experimental data and then proceed with 
the calculations within the dilution-frustration model of 
Eq. (10). Comparison with experimental data provides 
a confirmation of the selfconsistency of our considera- 
tion in two respects. First, it allows us to verify whether 
the impurity-induced frustration is at all a reasonable 
mechanism of producing extra suppression of the order 
parameter. Second, we are able to compare the range of 
parameters J' Zn and J^ n that are necessary to explain the 
experimental data for M(x) with the range permitted by 
the three-band Hubbard mapping of Sec. [Hi 



A. Discrepancy with experiments 



Our Fig. 10 'a) shows a comparison of the results for 
the averaged on-site magnetic moment (staggered mag- 
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FIG. 10. (Color online) (a) The on-site magnetic moment 
M(x) per magnetic site (Cu 2+ ), normalized to the undoped 
value M(0) vs Zn doping x. Green diamonds and magenta 
circles are the NQR datap^ cyan squares are the neutron 
scattering data™ and the dashed line is their best fit. The 
dotted line is the classical (Ising) result. Upper solid line 
(blue online) is the best fit of the QMC results (blue circles) 35 
and the lower solid line (red online) is the T-matrix calcu- 
lation resultsp' both for the dilution-only model, (b) The 
slope function R(x) of normalized staggered magnetization 
(see text) vs x. Symbols and lines are the same as in (a). 



netization) M(x), normalized to its value in the un- 
dope d sys tem M(0). The experimental data include the 
NQR^E^ and the neutron scattering datapH with the 
dashed line being their best fit. Both set of theoretical re- 
sults, from the QMCP^ and the T-matrhP^ calculations, 
are for T = dilution-only model, i.e., neglecting frus- 
trating effects of impurities proposed in this work. First 
observation is that the unbiased QMC data agree very 
closely with the T-matrix results up to x ~ 15%, thus 
supporting the validity of the latter in the low-doping 
regime pklMJ At higher doping, single-impurity T-matrix 
tends to overestimate the effect of impurities on the order 
parameter. Note that the T-matrix results in Fig. [lQ^a) 
are multiplied by the classical probability (dotted line). 
This does not affect the data for x < 35% but makes a 
comparison consistent near the percolation threshold x p . 

However, there are substantial discrepancies between 
theoretical and experimental results. In Fig. [lO^a) the 
experimental data are always below the theoretical curves 
for the dilution-only model. The difference is much more 
apparent in the slope function, defined as 



R(x) 



1 



1 M(x)\ 
x\ M(0) J ' 



(11) 



and displayed in Fig. [lQ^b). The quantity R(x) has a 
transparent physical meaning at x —> 0: it shows a re- 
duction of the order parameter M per impurity due to 
enhanced quantum fluctuations in the impurity's neigh- 
borhood, the quantity that should be captured properly 
by the T-matrix approach. Fig. lOjb) shows a large 
discrepancy of the experimental slope, i? exp (0) ~ 1.1, 
with both the QMC, R QMC (0) « 0.5, a factor of ap- 



proximately 2, and the T-matrix, R T (0) ~ 0.7, a fac- 
tor of about 1.6. This gives a clear indication that the 
dilution-only theory significantly underestimates the im- 
pact of the individual impurities on the quantum spin 
background. Thus, the dilution-only model is not enough 
to describe I^Cui-^Zn^O^ 



B. Qualitative Ji ff -Jf ff -Jf ff mean- field consideration 
of the dilution-frustration model 

We would like to argue that the longer-ranged spin 
interactions in the undoped system are unlikely to provide 
a resolution to the observed discrepancy in M(x). On the 
other hand, the dilution-frustration model of Eq. (10) 



can be approximated on a mean-field level as an effective 
model with further exchanges that are proportional to 
doping. Such a mean-field consideration offers a simple 
way of checking the viability of this model in explaining 
the discrepancy. 

One of the natural ideas to explain the disagreement 
is by including the longer-range interactions, such as J2, 
J3, etc., and ring-exchange, in the low-energy spin-only 
model for the undoped Cu02 plane. Such terms are 
generally present in any realistic low-energy spin mod- 
els of Mott insulators and formally come as a result of 
the higher-order expansion of the Hubbard-like models 
Although for the cupratesSSl sucn terms are of the or- 
der of 10% of the nearest-neighbor J, they do lead to a 
reduction^ of M. Qualitatively, assuming the presence 
of J2 term in the spin-model and using an expansion of 
M in the dilution fraction x and in J2, one can obtain 
on the mean- field level at small x and J2: 



R(J 2 )~R(0) [1 + A 



J 



(12) 



where the coefficient A<1 follows from the consideration 
of the J1-J2 modeP^ and R(0) is the theoretical slope of 
the dilution-only model in Fig. [l2^b). Thus, unless J2 
is of the same order as J, one can not expect a large 
correction to the slope of M(x)/M(0) from the presence 
of further-neighbor terms in the low-energy model. 

In addition, since these terms suppress the order al- 
ready in the undoped system, this mechanism is unlikely 
to enhance fluctuations specifically due to dilution. Put 
another way, dilution of the models with extended in- 
teractions breaks frustrated bonds alongside the nearest- 
neighbor ones and may even induce more robust order 
around impurities. A recent detailed study 38 has shown 
that while extended interactions such as J2 and cyclic 
terms do lead to an overall lower absolute values of the on- 
site magnetization in the impurity-doped Hubbard model 
of I^Cui-^Zn^O^ they do not significantly change the 
theoretical M(x)/M(0) dependence and are not able to 
explain the large experimental slope R(x) in Fig. 12 [b). 



On the other hand, the dilution-frustration model of 



Eq. (10) offers an alternative: while the undoped system 
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can be considered as having no further-neighbor frus- 
trating terms in its low-energy model, such terms are 
introduced by the dopants. We thus propose an "effec- 
tive medium", mean-field Jf ff -J| ff -J| ff model, in which 
effective interactions depend on the doping concentra- 
tion x, as a simplified version of the dilution-frustration 
model ( [To] ) that allows for a straightforward calculation 
of M(x)without the complications of the T-matrix or 
QMC numerical approaches. While approximate, this 
method gives an intuitive picture of the proposed mecha- 
nism and allows us to estimate whether it can be a viable 
source of enhancing quantum fluctuations due to doping. 

In a sense, the mean- field approximation "spreads" 
the frustration provided by impurities evenly over the 
whole system. The strength of the effective couplings 
is related to the exchanges in (10) by counting bonds: 
nearest-neighbor interaction jf^^= J(l — 2x), next- 
nearest-neighbor interaction J| ff = 2xJ' Zn , and next- 



Teff 



xJ z ' n . Since 



next-nearest-neighbor interaction J 3 
we are interested in x <C 1 limit, the derivation of M 
amounts to on expansion of the spin-wave correction to 
the order parameter within the Jf ff - J| ff - J| ff model to the 
first power in J| ff / J and J| ff /J. Using the approach of 
Ref . [56] we obtain the moment reduction as 



M(x) 
M(0) 



C 2 



c 3 



M(0) V J J Af(0) V J J ' 



Teff 
J 2 



(13) 



where C 2(3 ) = Ek ^2(3),k7k/ 2cj L with = a/T 

7k = (cos/c^ + COs/^)/2, 7^ 2 ,k = (1 ~ COS k x COS fc^) , 

^3,k = (2 - cos 2 fea. - cos 2 k y ), and M(0) = 0.3034 is the 
ordered moment of the 5 = 1/2, square-lattice nearest- 
neighbor Heisenberg antiferromagnet in the spin-wave 
approximation. Interestingly, since the C 2 = 0.2909 and 
C 3 = 0.5205, the suppression of the order by J| ff is al- 
most two times as effective as by J| ff of the same value. 
Setting the change of the slope, which is needed for the 
dilution-only result R T (0) « 0.7 to achieve the agreement 
with experimental i? exp (0) ~ 1.1, and using the relation 
of J| ff and J| ff to J' Zn and J Zn we can determine the 
range of J|°* = 27^+47^ necessary for that. Assuming 
J| ff = J| ff corresponds to the choice J Zn — 2J Zn , which 
is a reasonable ratio according to the three-band model 
consideration, see Fig. [8] Using (13), the required total 
frustrating effect due to each impurity is J z °^ ~ 0.6 J, 
which is well within the range of the estimates from the 
three-band model mapping, see Fig. [9] 

Altogether, this consideration shows that the dilution- 
frustration model is indeed a viable candidate for explain- 
ing the discrepancy with the experimental data. 



C. T-matrix for the dilution-frustration model 

A more rigorous, albeit more technically involved 
method of dealing with the dilution-frustration model 
( 10 ) is the T-matrix approach. Its advantage is in taking 



into account all multiple-scattering magnon processes, 
by which it solves the single-impurity problem exactly 



within the 1/5 spin- wave approximation. Since we are 
interested in identifying a mechanism of the order pa- 
rameter suppression in the low-doping regime, such an 
approach should be able to provide if not the ultimate, 
but at least a quantitatively correct result. 

The T-matrix method has been used in the study 
of the dilution- only model in the 2D square-lattice 
antiferromagnet pSHSU jj ere we outline some of the ba- 
sic steps needed to calculate staggered magnetization in 
the dilution-frustration model ( 10 ) and provide necessary 
details in Appendix [Bj 

Staggered magnetization at T = in the presence of 
impurities can be written as^ 



M{x) = M(0) 



E^(<4«k>-7k<4al k >) ,(14) 



where M(0) = S — A ~ 0.3034 is the staggered magne- 
tization of the 5 = 1/2, square-lattice, nearest-neighbor 
Heisenberg antiferromagnet at zero-doping, which is al- 
ready reduced by quantum fluctuations (A) of the anti- 
ferromagnetic ground state, and uj\^ = y/l — 7^ as be- 
fore. Quantum fluctuations due to impurities further 
reduce M(x) by non-zero magnon expectation values 
( a 3c a k) an d ( a ic a -k)- These can be expressed through 
the magnon Green's function, 34 leading to 



M(x) 
M(0) 



M(0) 



' duo 

00 f w k 



(imGiV) 



- ytfmG?{w)) , (15) 



where and G^ 2 are the diagonal and the off-diagonal 
component of the 2x2 matrix Green's function, respec- 
tively, see Eq. (Bl) for the corresponding Dyson-Belyaev 



expression of them via the magnon self-energies. We note 
that lmG~^(uj) > for uj < 0. Magnon self-energy in this 
approach comes as a result of averaging over random 
impurity distribution, which relates it to the forward- 
scattering (k 7 = k) elements of the T-matrix 



(16) 



Here the sum includes the contributions from impurities 
in both sublattices (I = A, B) and from each of the 
non-zero \i = s-, p-, and d-wave components of the T- 
matrix for the vacancy in the nearest-neighbor square 
lattice. Clearly, the self-energies are also proportional to 
the impurity concentration x. 



The individual components of the T-matrix (B5) obey 



linear integral equations for the multiple-scattering in 
each partial wave. Such equations contain separable scat- 
tering potentials, which reduce integral equations to sim- 
ple algebraic ones and allow for solutions in a separable 
form. That is, each partial wave T-matrix T^ /i k;/ (cj) is 
given by the product of the k- and k'-dependent scat- 
tering potential components with the c orre sp ondi ng in- 
dependent resolvents, 34 see Eqs. (B5), (B6), (B7), and 
jB8|) in Appendix pi 
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The dilution-frustration model (10) contains extra 
frustrating antiferromagnetic terms compared to the 



dilution-only model, J Zn 



and J Zn , that couple spins 



around impurities. We find that because these extra 
terms connect spins that belong to the same antiffer- 
omagnetic sublattice of the host, additional scattering 
provided by them is orthogonal to the s-wave due to a 
cancellation between S Z S Z and S+S~ contributions. For 
the same reason, the next-next-nearest-neighbor inter- 
action J Zn also does not modify the d-wave component 
of the scattering potential and contributes only to the 
p-wave, while the next-nearest-neighbor J Zn interaction 
affec ts bo th t he <i-wave and the p-wave scattering, see 
dB6l) and (|B8). 



While the x-dependence of Eq. (15) obviously extends 



beyond linear, it is the linear term which is of primary 
importance as it is defining the initial slope of the nor- 
malized magnetization, R(0) in (11). It is also the term 
that is expected to be treated properly by the T-matrix 
approach. Thus, one may want to simplify the general ex- 
pression (15) and obtain an expression which is explicitly 



linear in x. Having in mind that all magnon self-energies 
in (16) are ex x, we use the Green's function expansion 



in E, see (B2), to obtain 



M{x) 
M(0) 



1 - xR(0) = 1 



M(0) 



E 



-TkReS^(afc) 
2u£ 



J dw 

oo 



ImS^fw) 



7kIm£i 2 H 



) 2 (> 



,(17) 



where we note again that ImEj c 1 (cj) > for u < 0. With 
this, we can study the effect of impurity-induced frustra- 
tion on the order parameter by calculating integrals in 
(17) and (15) for different values of J Zn and J Zn . 

First, in agreement with the qualitative considera- 
tion provided by the mean- field J® ff -J| ff -J| ff model in 
Sec. Ill B, we find that the next-next-nearest neighbor 
J Zn bond suppresses the order about as effectively as two 
next-nearest J Zn bonds of the same strength. As was dis- 
cussed in our previous workp^l the same result has also 
been obtained by the QMC calculations, confirming again 
a good qualitative and quantitative accord between the 
analytical T-matrix approach and an unbiased numeri- 
cal method. Since the three-band model calculations of 
Sec. [Tl] also suggest that the J Zn -term is systematically 
larger than the J Zn -term, these findings make the J Zn 
interaction particularly important. 

We study the rate of suppression of the order parame- 
ter R(0) in (17) for several representative ratios between 

TT 



and J Zn as a function of one of them. Our Fig. 



shows R(0) vs J Zn for three choices of Jzn/^Zn = 0, 
and 1/2. The dilution-only T-matrix result 34 of R T (0) « 
0.69 and the "target" experimental value R exp (0) « 1.1 
are shown by the dotted horizontal lines. We find the 
effect of frustrating interactions to be about two times 
stronger than in the mean-field Jf ff -J| ff -J| ff considera- 
tion of Sec. Ill B, implying that the averaged mean- field 
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O Cu 
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FIG. 11. (Color online) R(0) from (17| vs J z ' n (in units of J). 
The red solid, magenta dotted-dashed, and blue dashed lines 
are for Jzn/^Zn — 0? 1/4 and 1/2, respectively. Colored area 
highlights the range of J z ' n . The experimental value R exp (0) 
and the dilution-only T-matrix result 34 of R T (0) are shown 
by the dotted horizontal lines. 



treatment underestimates the effect of the order suppres- 
sion due to frustration by local defects. This result also 
means that the necessary frustration is smaller than the 
one estimated from the Jf ff - J| ff - J| ff model. 

One can see that the experimental value of R(0) in 
Zn-doped La2Cu04 is met by the T-matrix results of the 
dilution-frustration model at J Zn « 0.07 J if the ratio of 
J' Z J J z ' n is fixed to 1/2 and at J z ' n « 0.08 J for J Z J J z ' n = 
1/4, the window of variation of T Zn /J Zn suggested by 
the three-band model consideration in Fig. [9] of Sec. [TTJ 
These values translate into the total frustrating effect 
of impurity J z ° £ w 0.28 J and J|°* w 0.24J, respectively, 
which are well within the window suggested by the three- 
band model calculations. The somewhat wider range is 
highlighted by the gray shaded area in Fig.^a) and (b). 

In our previous workp^ QMC results seem to sug- 
gest a higher value J|°* > 0.4 J, still a modest amount 
of frustration, well within the range permitted by the 
three-band model. An example of the QMC results for 
R® MC (0) for the choice of parameters J z ' n = 2J Zn = 0.1J 
is shown in Fig. [l2^b) by the blue diamond on the ver- 
tical axis. However, the QMC calculations may be af- 
fected by the frustrating nature of the impurity-induced 
interactions, which is associated with the infamous sign 
problem. This is also restricting the use of the QMC 
for the dilution-frustration model to the single-impurity 
problem. In addition, the QMC results for R® MC (0) are 
obtained from the finite-size scaling that may overesti- 
mate Jrjff . Thus, using the data only from the largest 
clusters,^ extrapolation for the J^ n — 2J' Zn — 0.1 J data 
set gives R® MC (0) « 1.0, much closer to the experimental 
value. Thus, the QMC result also demonstrates that the 
impurity-induced frustrations affect the staggered mag- 
netization significantly and bring the slope much closer 
to experimental data. 
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QMC, dilution-only 
J-matrix, dilution-only 
r-matrix, 7 Z '„'= 2/ z '„ =0.07 
QMC, J^= 2 j; n =0.1 
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FIG. 12. (Color online) Same as Fig. [To] with the lines used 
for the fit of the QMC data and for the T-matrix results of 
the dilution-only case switched to the dashed lines for clar- 
ity. The red soli d lin e (and its long-dashed line tail) is the 
T-matrix results (15) for the dilution- frustration model with 
Jz n = 2 = 0.07J. In (b) the blue diamond with the dashed 
line is the QMC result for j£ n = 2J' Zn = 0.1J from Ref . l39l 



Using Eq. (15) we calculate the normalized stag- 
gered magnetization and the slope ( pTj ) for the dilution- 
frustration model with J% n = 2J! Ln = 0.07 J as a function 
of with the results shown in Fig. 12 'a) and (b) by the 



solid lines. The comparison is provided with the experi- 
mental data and the theoretical T-matrix and QMC cal- 
culations for the dilution-only model from Fig. [10] With 
the frustrating parameters chosen to match the initial 
slope of the experimental best fit, obviously, M(x)/M(0) 
in the dilution-frustration model agrees much better with 
experiments. At higher values of x > 15%, the T-matrix 
overestimates the effect of impurities on M(x) as it ne- 
glects the multiple-impurity scattering effects, similarly 
to the dilution-only case. 

Altogether, all of the evidences indicate that the 
dilution-frustration model describes La2Cui_ :c Zn :E 04 
better than the dilution-only model and that the effective 
frustrating interactions via available electronic states of 
the impurity should be taken into account when studying 
doped Mott- and charge-transfer insulators. 



D. Further developments 

Based on the analysis of the Neel temperature dop- 
ing dependence within the classical Heisenberg model, 
an alternative suggestion has been put forward for ad- 
ditional interactions introduced by dopants, such as Mg 
and Zn: lattice distortions in the vicinity of the impu- 
rity site change the strength of J bonds.^The impact of 
this mechanism on M(x) and its possible viability as an 
alternative to our proposal has been investigated using 
the QMC analysis in Ref. [58) It was found that changing 
the strength of the eight Cu-Cu bonds in the immediate 
vicinity of Zn site by as much as 15% (= SJ to t — 1.2 J) 



changes the slope of M(x)/M(0) by at most a few per- 
cent. Such a weak effect rules out the lattice-distortion 
mechanism as a viable alternative to our theory. 

One of the natural consequences of our idea is that 
different isovalent impurities, such as Zn 2+ and Mg 2+ , 
should induce different amount of quantum fluctuations 
around them due to differences in their electronic levels 
available to generate frustrating interactions, and thus 
lead to a different rate of suppression of the order pa- 
rameter. Specifically, Mg 2+ should have no levels in any 
reasonable proximity to the chemical potential, as op- 
posed to Zn 2+ , and thus should not lead to any substan- 
tial frustration, suggesting that the Mg-doped La2Cu04 
must be described much better by a simple dilution-only 
model. This scenario has been investigated recently by 
the /iSR experiments in Zn- and Mg-doped La2Cu04 at 
low doping™ The crucial finding of this work is that 
the spin stiffness, determined from the Neel tempera- 
ture, shows a stronger suppression rate in the case of 
Zn-doping, in agreement with the expectations from our 
theory. Although in additional to quantum contributions 
the slope of T/v(x)/T/v(0) « 1 — ax contains significant 
classical terms as well as logarithmic contributions from 
the three-dimensional interplane couplings j 34 l 57 l the dif- 
ference in such slopes found in Ref. 40 is significant: 
«M g ~ 2.7 vs azn ~ 3.5. These values are consistent 
with the difference expected between the dilution-only 
model 34 and the dilution-frustration model with the pa- 
rameters discussed in Sec. IIIII C. 



IV. CONCLUSIONS 

In this work, we have provided a detailed deriva- 
tion of the low-energy, spin-only model of the Zn-doped 
La2Cu04 starting from the realistic site-diluted three- 
band Hubbard model. We have followed with an equally 
detailed analysis of the order parameter in this low- 
energy model. We have elaborated on the proposal of 
our previous work that the impurities in strongly cor- 
related systems may induce significant longer-range frus- 
trating interactions among the spins, which are absent or 
negligible in the corresponding undoped system, due to 
hybridized electronic states of the impurity at the scale 
less than the Hubbard- U. Such impurity-induced frus- 
trating interactions are shown to significantly enhance 
local quantum fluctuations compared to the dilution- 
only model. In particular, the dilution-frustration model 
demonstrates stronger suppression of the order and, for a 
choice of parameters appropriate for Zn-doped La2Cu04, 
it resolves discrepancies between experiments and ear- 
lier theories. Recently, our theory has received further 
experimental confirmation from the /iSR studies of spin 
stiffness in Zn and Mg doped La2Cu04. 

Although our work considers a particular case of di- 
luted La2Cu04, this study has far-reaching consequences 
for other diluted antiferromagnets and doped Mott and 
charge- transfer insulators One of the intriguing conse- 
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quences of the proposed mechanism is the change of the 
character of the percolation transition due to frustrating 
interactions across the impurities. The impurity doping 
of one-dimensional spin systems such as spin chains and 
ladders should introduce weaker links between their parts 
instead of breaking them into independent pieces. Exper- 
iments in diluted frustrated spin systems is yet an other 
area, such as recently studied diluted J1-J2 system ) 60 * 61 ! 
which should also be affected by the same mechanism. 
Another perspective is offered by an extension of the pro- 
posed mechanism of the impurity-induced frustrating in- 
teractions to the case of the doping away from half-filling 
where it may be responsible for the pair-breaking mech- 
anism in the doped Cu02 planes. 
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the Cu02 plane doped with Zn, the Fourier transform of 
the oxygen-hole operators p™^ can be written as 



1 rncx 



l=X>- k -(n ± f(¥)) p g), (A3) 



where the sum is over the Brillouin zone and N is the 
number of unit cells. 



Motivated by (A2) and using (A3), one can introduce 



the orthogonalized symmetric and antisymmetric oxygen 
operators in momentum space as 



cos I y^Pk,a + cos b } L 

COS( yj a - COS ) PkA 



(A4) 



with the normalization factor 



■7k, 



(A5) 



and 7k = \ (cosk x +cosk y ). These operators are now 
properly normalized and obey regular fermionic anticom- 
mutation relations. 
Then, using 



k 



-ik r; \ 
A 



k #k,o 



(A6) 



and applying the inverse Fourier transform yields the 
hopping Hamiltonians in Q and ([5| with the real-space 
Wannier amplitudes given by 



ke 



ik-(rz-ry) 



(A7) 



Appendix A: Details of the three-band model 
consideration 

1. Wannier orthogonalization of oxygen orbitals 

The hopping terms in the three-band model and ([2| 
couple the Cu and Zn states with the symmetric combi- 
nation of O-states, i.e., 



with 



H pd = -t pd Y,(dlP la +K.c. 



(Al) 



■) , (A2) 



where x and y are the vectors in the x and y di- 
rection of length ao = 1, the Cu02 lattice spacing. 
The problem is that the operators Pi a do not obey 
proper anticommutation relations with the ones on the 
nearest-neighbor CuO ^Zn O^ clusters, because one of 
the four p-operators in ( A2 ) belongs to both clusters. The 



Wannier-orthogonalization procedure is performed in the 
k- space!- 43 * 46 * Since the O-lattice is fully periodic even for 



which depend on the distance between the clusters. As 
a result, the original three-band Hamiltonian ([TJ with 
Zn-impurity Q is rewritten in a much more convenient, 
symmetric CuO^ZnO^ cluster form, containing oxygen 
degrees of freedom that couple to copper and zinc orbitals 
most effectively. On the other hand, the hopping terms 
in the new variables go beyond the nearest-neighbor hop- 
ping. However, \n> decrease rapidly with distance,^ 
\n> ~ 1/1 — r/' | 3 , so the terms beyond the nearest neigh- 
bor in the hopping part of Q and ([5| can be neglected. 

The most important result of the transformation 
{Pm,ai Pma} ~^ ^,afe,a) 5 is that it allows to take into 
account the effects of the the intra-cluster hoppings sep- 
arately from the local, intra-cluster Cu-0 (Zn-O) hy- 
bridization (Wannier amplitude Aq), considered next. 



2. Effective CUO4 and Zn04 states and hoppings 

Here we list a complete set of hybridized orthogonal 
electronic states in CUO4 and Zn04 clusters for the one- 
and two-hole states. We omit the site-index as the states 
are in the same cluster. The states involving antisym- 
metric oxygen orbitals (A4) are excluded since they do 
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not contribute to the hybridization. With zero-hole state 
defined as |0), the CuC^ one- hole states include 



K) = 4|0), \q a )=qt\0), 



(A8) 



with energies Sd and £ p , respectively. 

The two-hole sector is naturally divided into orthogo- 
nal singlet and triplet sectors, to be diagonalized sepa- 
rately, with singlet states 



|V)=44|0), \cp)=q\ql\0) 

lx) = 71(4^-44)1°)' 



(A9) 



having energies 2s d + Ud, 2e p , and Sd + £ P , respectively, 
and triplets 



\r 1 ) = d\q\\0), |t_i)=4 9 J|0>, 
1 

Ti' 



ir > = -W44 +44)io), 



(A10) 



all with the energy Sd + s p . 

Similarly, for ZnC>4 cluster the one-hole states are 



K) = 4|o), |« a > = «i|o), 



(All) 



with energies sz n and e p , respectively. The two- hole sin- 
glet and triplet states are 



iv Zn > = 44io), 



\^ n )=q\q\\0), 



ix z -) = ^=(44 -44)io), 

with energies 2e Zn + /7 Zn , 2£ p , £ Zn + £p, and 



(A12) 



|rD=44|0), \T^) = a[ql\0) 



t„ti 



10 = ^(44+44)1°)' 



(A13) 



with energy ez n + s p . The states involving only oxygen 
orbitals are the same for both Zn04 and CUO4 clusters. 

Hybridization terms in 7-L loc in Q and in m @ 
originate from hoppings between the orbitals within each 
cluster and lead to mixing of th e st a tes in each of th e 
ortho gonal sectors of states in (A8)-(A10) and (All)- 
(A13). Thus, the diagonalization of the local Hilbert 



space in each cluster amounts to bringing a few 2x2 
and 3x3 matrices to a diagonal form. We do not list 
explicit expressions for the resulting eigenenergies and 
eigenvectors, which can be found in Ref. [46l and simply 
assume that they can be easily determined for a given 
choice of the three-band model parameters. 

Since we are interested in the lowest states from each 
of the n-hole sector, we denote them as follows. The 
lowest one-hole state of the Cu-cluster is \fa^), which 
is a normalized linear combination of \d a ) and \q a ) in 
(A8), with the energy E\. The lowest one- hole state for 
Zn- clu ster is 1/^), a combination of \a a ) and \q a ) from 
(All), with the energy E\. The lowest two- hole states 
are \ f^) for Cu-cluster, which is a combination of the 



singlets \ip), and \x) in (A9), with the en ergy E 2 and 
|/v 2 )) f or Zn-cluster, a mix of the singlets in (A12), with 
the energy E 2 . 

Setting the energy E\ of l/^^-state on Cu-cluster to 
zero, as in Figs. [3] and [4j defines the relevant energy scales 
of the effective t-e-U model as = E\ — E\ for the 
lowest one- hole state of Zn-cluster, U^ 1 = E 2 —2Ei as the 
effective Hubbard gap on Cu-cluster, and = E 2 -2E 1 
for the effective repulsion on Zn-cluster. Ignoring the 
hoppings that involve higher-energy states, the effective 
t-e-U model, which is an abbreviated version of the model 
in ([6| and (|7|), is given by 

^ = E C/ e C ffl// 2) )(// 2) | (AM) 

I 

E{^ 1 (lil! ) )|o i ')e ) l(o 1 | + H.c.) 
(II') 

F$ (|/^)|/^)(0H<// 2) | + H.c. + « I'}) }, 
and 

= E { m + ^eDl/PX/Pl + E^l/^)^!} 
i a 

+ E{ i? i°o (l/^ ) )|O i )(/,il ) KO.|+H.c.) (A15) 



(a) 



+ ^?(l/^ ) >l/iJ ) ><0iK/7 ) l+H.c.) 
+ ^i 2 (l/^)l/^)<// 2) K0.|+H.c.)}, 



where the explicit expressions for hopping integrals can 
be obtained by evaluating matrix elements of the Hamil- 
tonian in terms of the original Cu, Zn, and O operators, 
Q and ([5]), between the initial and final states in the 
basis of the local eigenstates. For example, 



77101 
^10 



= (/S ) I(OH(w + *«)|/^ ) >|Oz>- (A16) 

In the effective hopping terms shown in Fig. [5] and used 
in Sec. JIT] we have employed the shorthand notations: 
^1 = ^10 > ^2 = F u , tn = F 10 •> ^21 = ^11 •> and ti2 = F n . 



Appendix B: Details of the T-matrix calculation 

Here we provide some of the technical details of the 
T-matrix approach, which follows closely Ref. [34j 

The 2x2 matrix Green's function of a magnon in the 
square-lattice antiferromagnet can be written in a Dyson- 
Belyaev form: 



Gk(w) 



-1 



-oj — CJk - 



^22 



UJ — CJk — 



(Bl) 
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The low-doping consideration requires an expansion of 
(Bl) in powers of E, which given by 



G ll = G 0,11 + G 0,11 E 11 G H + G 0,11 E 12 G 21 

- G ' 11 + G ' 1 ^ 11 ^ ' 11 + 0(x 2 ) , (B2) 

q12 _ qQ, 11^11^12 _p gfO, 11^12^22 

^G°' n E 12 G°' 22 + (9(x 2 ) 

where we drop the common k and u dependencies for 
shorthand notations and the non-interacting magnon 
Green's function is 



i 



V 



■ w k ■ 
o 



(B3) 



U + Co>k — zOy 



where we normalize all the energies to = 45J and 
the magnon frequency is uj^ = y^l — 7^, with 7k = 

(COS fc^ + COS fc^) /2. 



The self-energies in (Bl) and (B2) are linear in the 



doping concentration x and are related to the forward- 
scattering components of the T-matrix via 



(B4) 



where the sum is over the sublattice index (7 = A, 5) 
and /i = s-, p-, and d-wave components. 

Using the sublattice-A as an example,^ the s-, p-, and 
d-wave components of the T-matrix are 

^H = l) ^r,H (M€p,d), (B5) 
+I^)0(a4i + |a^)0(4i, 

where we use <g> to denote the direct product of the col- 
umn and row vectors and V^w stand for the s-, p-, and 
d-wave components of the scattering potential that can 
be written as 



fj,Es,p,d 



(B6) 



in which the s-, p-, and d-wave vectors are given by 

(«k I = (w k , -Vk) , (As£| = K, Vk) , (B7) 
(Pitc) I = ^ sin ^(j/) K, Uk) , (dk I = 7k K, ^k) , 



with the Bogolyubov parameters ^k = \/(^+ ^k)/2cJk 
and ^k = — sgn(7k) \/ (1 — Uk)/2^k- The cj-dependent 
resolvents r /i (cj) in ( |B5|) are listed below, see (B8). 

The coefficients in ( |B6| ) contain the dependence on the 
frustrating terms in the dilution-frustration model 



G s = l, G p = l + 2^+2^, C a 

For the sublattice- 7?, the equivalent expressions are ob- 
iined viaf^(cj) = 
The resolvents are 

(1 + lj) P (lj) 



v 

1 + 4^ 



tained via T^£,(uj) = T^,(—uj){u <R> i?}, see Ref. 



r s (w) 
r p (w) 
r rf (w) 



w(l + w)p(w) 



2 - CP (1 - w) [1 - w 2 p(w) + Pd (u)] 
-1 



(B8) 



l + c^i-^pdH ' 

where the functions p(u>) and Prf(w) are 



m^) = E 



(B9) 



with 7^ = (cosk x — cos k y )/ 2. They can be expressed 
through the complete elliptic integrals of the first and 
second kindP^ 

The self-energy matrix elements are related via 
E 22 (cj) = E^(-cj) and E^ 2 (-cj) = E 21 (u;). The partial- 
wave terms in the self-energies are 



xC s uj h 3 v 7 V^ k 1 



for the s-wave, 



cj k 




■r-H 



for the p-wave, and 



xC d w k 



2k 



'{rj(") 



^ V° - 1 
1 
1 



1 -7k 

-7k 1 

-o; k 
o; k 



1 -7k 
-7k 1 



(BIO) 



}■ 



(Bll) 



-w k 



J}' (B12) 



for the d-wave. Here = - [^(cj) ± cj)] and 
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